function [beta,tOls,tWht,tNwest,r2,adjr2]=ols(X,Y,lag)
[N,k] = size(X);
beta = X\Y;
e = Y-X*beta;
ebar = Y-mean(Y);
Q = inv(X'*X);
r2 = 1 - (e'*e)/(ebar'*ebar);
adjr2 = (1-k)/(N-k) + (N-1)*r2/(N-k);
varOls = ((e'*e)/(N-k))*Q;
tOls = beta./sqrt(diag(varOls));

varWht = e.^2;
varWht = X'*(X.*varWht(:,ones(k,1)));
% HTR consistent variance = inv(X'X) (X'SX) inv(X'X)
%   multiplication by N/(N-k) improves small sample properties
varWht = Q*varWht*Q*N/(N-k);
tWht = beta./sqrt(diag(varWht));

% Newey West consistent variance ref Campbell,Lo,Mackinlay pp. 534
%lag = floor(4*((N/100)^(2/9)));
v = e.^2;
v = X'*(X.*v(:,ones(k,1)));
for i=1:lag
    v1 = e(1+i:N).*e(1:N-i);
    v1 = X(1+i:N,:)'*(X(1:N-i,:).*v1(:,ones(k,1)));
    v = (1-i/(lag+1))*(v1+v1') + v;
end
varNwest = Q*v*Q*N/(N-k);
tNwest = beta./sqrt(diag(varNwest));
return